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We study the non-equilibrium steady state arising from the interplay between coherent and dis¬ 
sipative dynamics in strongly interacting Rydberg gases using a recently introduced variational 
method [H. Weimer, Phys. Rev. Lett. 114, 040402 (2015)]. We give a detailed discussion of 
the properties of this novel approach, and we provide a comparison with methods related to the 
Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy. We find that the variational approach offers 
some intrinsic advantages, and we also show that it is able to explain the experimental results 
obtained in an ultracold Rydberg gas on an unprecedented quantitative level. 
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I. INTRODUCTION 

The discovery of dissipative quantum state engineer¬ 
ing [InSl, i' e ’’ th e use of controlled sources of dissipation 
for the preparation of many-particle quantum states, has 
led to a surge of interest in open quantum many-body 
systems. The unrivaled tunability of interaction and dis¬ 
sipation properties of driven Rydberg gases makes them 
particularly useful for this purpose and has resulted in 
a large amount of theoretical and experimental works 
investigating the interplay between coherent and dissi¬ 
pative dynamics in these systems UH3]- Here, we pro¬ 
vide a detailed discussion of the application of a recently 
introduced variational principle for the steady state of 
dissipative quantum many-body systems [25} to driven- 
dissipative Rydberg gases. 

Rydberg atoms are routinely excited from ground state 
atoms by coherent laser driving, with their atomic prop¬ 
erties and interaction strength scaling dramatically with 
the principal quantum number |26>]. The radiative decay 
of the metastable Rydberg state provides a natural dissi¬ 
pative element whose decay rate can be widely tuned by 
laser coupling to other excited non-Rydberg states [271 ] . 
As such, Rydberg atoms provide an ideal environment for 
studying dissipative many-body dynamics in a strongly 
interacting regime. 

In this article, we perform an analysis of the non¬ 
equilibrium steady state of a driven-dissipative Rydberg 
gas. We investigate the properties of this steady state 
using a variational method recently introduced by the 
author [25[. We provide a detailed evaluation of the ap¬ 
proximations carried out within this novel method. We 
complement our analysis of the variational method by 
describing an alternative approach to the investigation 
of dissipative quantum many-body systems based on a 
hierarchy of equation that also allows for a systematic 
incorporation of correlations. We perform an explicit 
comparison of the two methods, finding substantial ad¬ 
vantages in favor of the variational approach. Finally, 
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we demonstrate that the variational method provides re¬ 
markable quantitative agreement with the experimental 
results on the steady state phase transition observed in 
an ultracold Rydberg gas. 


II. DRIVEN-DISSIPATIVE RYDBERG GASES 


We first give a microscopic description of driven- 
dissipative Rydberg gases. We express the dynamics in 
terms of a spin 1/2 representation, where the down spin 
state corresponds to an electronic ground state, and the 
up spin state refers to a highly excited Rydberg state [2(| . 
As the frequency difference between all electronic states 
is much larger than the decay rate of the Rydberg state, 
it is well justified to describe the radiative decay of the 
Rydberg excitations in terms of a Markovian quantum 
master equation in Lindblad form, 
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where p is the density operator describing the state of the 
system, H is the Hamiltonian accounting for the coher¬ 
ent part of the dynamics, and the set of Ci are quantum 
jump operators responsible for the dissipation [281 ]. If the 
external laser fields are close to a resonance between the 
atomic ground state and a single Rydberg state, such a 
system in its rotating frame is governed by the spin 1/2 
Hamiltonian 
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where the Rabi frequency f 1 and the detuning A repre¬ 
sent the laser parameters, the Cq coefficient denotes the 
strength of the van der Waals repulsion between Rydberg 
states, and P[ = (1 + of )/2 is the projection onto the 
Rydberg state. The jump operators describing the decay 
of the Rydberg excitations is represented by quantum 
jump operators describing up spins flipping into down 
spins according to c, = y/jcr^}, with 7 being the decay 
rate of the Rydberg state. Initially, we will assume that 
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the system is well described by a model involving only 
nearest-neighbor interactions, i.e., the blockade radius is 
smaller than the lattice spacing a. Then, we can intro¬ 
duce an interaction constant given by V = Ce/a 6 . Of 
special interest is the regime where the Hamiltonian is 
equivalent to a transverse field Ising model, which is re¬ 
alized for A = V/2. 

In contrast to closed quantum systems, an open quan¬ 
tum systems such as the dissipative Rydberg gas will 
generically relax towards a stationary state character¬ 
ized by the condition 4rp = 0. Crucially, the inter¬ 
play between coherent and incoherent part of the dy¬ 
namics generically leads to a stationary state that is dif¬ 
ferent from any state in thermal equilibrium, i.e., a non¬ 
equilibrium steady state. 


III. VARIATIONAL PRINCIPLE 


exponentially with the system size, so additional steps 
are needed first. This situation is very similar to that of 
correlated fermions, where energy expectation values for 
variational Gutzwillcr wave functions can only be evalu¬ 
ated within further approximations 31]. In our case, we 
exploit the fact that we are not so much interested in the 
actual value of the trace norm, but rather in the proper¬ 
ties of the stationary state. The additional approxima¬ 
tions we will carry out retain the variational character of 
our calculation, i.e., they provide a rigorous upper bound 
to the trace norm. 


A. Product states 

To be explicit, we first consider the case of the vari¬ 
ational set of states being restricted to product states, 
i.e., 


In the following, we review the basic concepts behind 
the variational principle introduced in Ref. [251 ] . The ba¬ 
sic idea is to take a variational trial state p and compute 
the residual dynamics it will generate by computing its 
time-derivative p = Cp according to the underlying quan¬ 
tum master equation with the Liouvillian superoperator 
C. The operator p can be expressed as a traceless Hermi- 
tian matrix. To find the variational approximation to the 
true steady state having p = 0, we choose the variational 
state that minimizes the trace norm of p, i.e., 

Pvar = arg minTr{|£p|}. (3) 

p 

Choosing the trace norm as the correct matrix norm can 
be motivated on two different grounds. First, the trace 
norm is unbiased in the sense that it does not favor cer¬ 
tain classes of variational states over others without a 
physical reason. This is related to the linearity condition 
||p|| = ||Ap||/A satisfied by the trace norm. In particular, 
any other Schatten norm ||p|| = Tr{|p| p } with p > 1 is 
biased towards the maximally mixed state [25| . 

The second way to motivate the choice of the trace 
norm follows follows from quantum information theory. 
Here, the trace norm can be interpreted as being equiv¬ 
alent to the trace distance of p to the zero matrix, while 
the latter is obtained for p if and only if the variational 
state p is an exact stationary state of the master equa¬ 
tion. Importantly, the ability to physically distinguish 
the operator p from the zero matrix is given by their 
trace distance [29| . and hence the trace norm of p. In 
this sense, the trace norm is the natural norm to de¬ 
cide which variational state is the best approximation to 
the true steady state. We would also like to point out 
that the trace norm for the steady state is equivalent to 
applying a time-dependent variational principle for the 
dynamical evolution [30]. 

Following these initial statements, we can now pro¬ 
ceed with the variational analysis. Calculating the trace 
norm is in general still a computational problem scaling 


p = ni = Y[ Pi . ( 4 ) 
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Here, we have introduced the superoperator TZ, which re¬ 
places every occurrence of the identity operator for site 
i , 1 i, by the single-site density matrix pi. Addition¬ 
ally, we focus on quantum master equations including 
nearest-neighbor interactions or jump operators involv¬ 
ing at most two adjacent sites. Then, the trace norm of 
the resulting dynamics can be written in the form 

IIpII = IIE 7 ^ + E 7 ^H’ ( 5 ) 
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where pi = Tr ^{p} describes the single-site dynamics and 
Cij accounts for correlations between the sites. Here, 
the correlations between the sites stem from the nearest- 
neighbor interactions and two-site jump operators and 
thus are restricted to nearest-neighbor correlations only. 
As a first approximation, we apply the triangle inequality 
to pull the summation over i out of the norm, 

lldl <Eh^ + E^H- ( 6 ) 

* 3 

As the next step, we make use of the fact that pi and pj 
act on different parts of the Hilbert space, which allows 
us to write 

I Id I < E + E^ {pH + 3y) II- (7) 

* 3 

Note that this inequality does not change the variational 
approximation to the steady state, but it allows us to 
write the final result in a more compact form. Finally, 
we employ the triangle inequality a second time, yielding 

I Id I < E 11^ (dPj + PiPj +Cij) II = E Tr {idj|} ( 8 ) 

(u) <*J> 

Consequently, we have succeeded in mapping the full 
quantum many-body problem into a sum of efficiently 
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FIG. 1: Difference in the variational states introduced by the 
inequalities leading to Eq. (JHJ) as indicated by the trace dis¬ 
tance between the states p and p . Here the full space density 
matrix p was taken as a product state of four sites. The in¬ 
set shows the difference in the single particle reduced density 
matrices at = 4y (A = V/2, V = 5 7 ). 


solvable problems involving only neighboring sites. For 
translationally invariant systems it is in general even suf¬ 
ficient to solve a single two-site problem. 

At this point, it becomes a natural question to ask 
how well justified our approximations are. Besides the 
obvious approximation in restricting to the variational 
manifold, we have to assess the deviations introduced by 
the inequalities leading to Eq. (0. For small system sizes, 
we can answer this question exactly because the Hilbert 
space is still small enough so that we can minimize the 
trace norm of Eq. 0 without any additional approxima¬ 
tions. As shown in Fig. [l] the deviation introduced by 
the inequalities is quite small for the Ising model describ¬ 
ing the dissipative Rydberg gas and in the case of four 
particles. 

Of course, we are not really interested in problems in¬ 
volving four sites, but rather see how the situation be¬ 
haves in the thermodynamic limit. For this, we can check 
how the exact variational norm behaves as a function of 
system size. As a single evaluation of the norm is compu¬ 
tationally less costly than performing a full minimization, 
we can go to somewhat larger system sizes. Remarkably, 
we find that the scaling with system size is completely 
independent of the model being investigated or the trial 
state for which the norm is evaluated, see Fig. [21 For 
a simple non-interacting toy model, this surprising fact 
can be understood as a consequence of the central limit 
theorem. This model consists of N purely dissipative 
two-level systems, whose jump operators are given by a 
dissipative spin-flip of the form c* = • v / 70 '“, as in the 
case of the dissipative Rydberg gas. As a trial state, we 
choose the maximally mixed state, p = 1/2^. Since the 
model is purely classical and non-interacting, we can give 
an analytical expression for the trace norm of the master 
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FIG. 2: Trace norm of the master equation depending on the 
system size N for different models and different product states 
normalized to the value at N = 4. The asymptotic behavior 
for the noninteracting decay model is shown as a solid line. 


equation, 



This expression can be evaluated efficiently even for large 
values of N and corresponds to the crosses in Fig. [2] To 
obtain the asymptotic behavior in the limit of large N, 
we replace the sum by an integral and use the central 
limit theorem to approximate the binomial coefficients 
by a Gaussian function, obtaining 



As shown in Fig. [21 this asymptotic behavior is already 
reached for quite small values of N, indicating that the 
central limit theorem can also be applied to the fully 
quantum case, which appears to be a natural consequence 
of the correlations in p being restricted to nearest neigh¬ 
bors. Therefore, we conclude that the approximations 
needed for an efficient calculation of the variational norm 
are well justified and enable to use the variational princi¬ 
ple as a powerful tool to compute steady state properties. 


B. Correlations 

The considerations made for product states can also 
be extended towards more generic classes of variational 
states including correlations. Here, we study the case 
where nearest-neighbor correlations are fully included, 
resulting in a variational state according to 

p=n*+£ nc l3 + J 2 KC^Ckj + .... (11) 

i ( ij) 

Additionally, we impose the constraint that all reduced 
density matrices pij = piPj + Cij are positive definite. 
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Using the same steps as for product states and observing 
that partial traces of correlations vanish, Tr i{Cij} = 0, 
we can again find an upper bound to the variational norm 
as 

m\< (12) 

(ijk) 

Here, the many-body problem reduces to a sum of three- 
site problems. It should not be surprising that the in¬ 
clusion of all two-site correlations leads to the minimiza¬ 
tion of a three-site problem as the interaction generically 
generates one higher order of correlations. Consequently, 
treating n-body correlations exactly requires to solve an 
n + 1-body minimization problem. 


by going up to higher terms in the hierarchy. However, 
the main drawback of the method remains, that it can¬ 
not be formulated in terms of a variational principle, i.e., 
the neglected higher order terms are uncontrolled. In the 
case of the hierarchy equations having multiple solutions, 
it is actually possible to combine it with the variational 
method. First, all solutions to the hierarchy equations 
are computed, which are then used as a variational class 
of states to find the solution that leads to a minimization 
of the variational norm. 


V. RESULTS 


IV. HIERARCHY EQUATION METHODS 


A. Lattice model 


An alternative way to analyze dissipative many- 
body dynamics is through a hierarchy of equations 
in terms of their correlations, in close analogy to 
the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy of 
classical physics [32j. For this, we express the density 
operator p in terms of their reduced density matrices, 
according to the generating functional 


T(a) =logTr<pJJ(l. 1 + a i ) \ , 


(13) 


where the operators on fornr an arbitrary operator ba¬ 
sis acting on lattice site i [33j. The first terms of the 
hierarchy are then given by 


Pi = Tr v {p} = 


dT 

dcti 


a—0 


Pij — Tr i/y{p} — piPj + 


d 2 T 
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(14) 

(15) 
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Pijk — Tr^ {p} — PiPjPk + CijPk H - CikPj ~\~ CjkPi 
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The same analysis can also be performed on the level of 
the quantum master equation to obtain effective equa¬ 
tions of motion for the reduced density matrices. Gener¬ 
ically, each equation of motion is coupled to the next 
higher equation of motion in the hierarchy. The usual 
strategy is then to truncate the hierarchy at some point 
by setting the contributions from all higher order deriva¬ 
tives of T to zero and solve the resulting closed set of 
equations 01- This approximation is attributed to a 
\/z suppression of the higher order derivatives. In low¬ 
est order, the resulting equations of motion are identical 
to the mean-field decoupling of Ref. [Hj. 

In principle, it is possible to systematically incorpo¬ 
rate correlations similar as in the variational approach 


We will first put our attention to the case where the 
atoms are distributed on a two-dimensional square lat¬ 
tice. We will make a direct comparison between the vari¬ 
ational method and the results from solving the hierarchy 
equations. We complement this comparison with results 
from a numerical solution of the quantum master equa¬ 
tion using a quantum trajectories method j35| . To ensure 
a comparison on an equal footing, we will compare the 
variational results for product states to the first order 
hierarchy equations, and the variational method for cor¬ 
related states to the second order hierarchy equations. 
In the latter case, we include only nearest-neighbor cor¬ 
relations within both methods. The results are shown 
in Fig. [3] Except for intermediate values of Q, the two 
methods agree very well, and for correlated states, the 
level of agreement is further improved and also matches 
well with the results from the quantum trajectories sim¬ 
ulation. 

However, some important qualitative differences re¬ 
main even when correlations are included. The bistabil¬ 
ity of the first order hierarchy solution is still present, 
although over a smaller range of parameters. Conse¬ 
quently, this bistability is not just an artifact of the first 
order result and the approximation of neglecting the 1 /z 
corrections to it. Rather, it appears to be a generic el¬ 
ement of the hierarchy equation method. In contrast, 
the variational solution always produces a unique sta¬ 
tionary state, as even in the case of multiple local minima 
of the variational norm, there is always a unique global 
minimum. Within the variational approach, we find a 
first-order phase transition between a low-density gas of 
Rydberg excitations and a high-density liquid 25:]. 

Finally, we turn to the parameter regime for nonzero 
detuning A where the mean-field solution (i.e, the first- 
order hierarchy equation) predicts the existence of an 
antiferromagnetic phase [5|. Within the variational 
method, we also find such an antiferromagnetic phase, 
see Fig. U but its extension is reduced significantly. 
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FIG. 3: Comparison of the Rydberg density n r for the variational solution (solid) and the solutions of the hierarchy equations 
(dashed) for product states (left) and correlated states (right). For reference, the quantum trajectories solution of the quantum 
master equation for a 4 x 4 system is shown as a dotted line. The data for both the variational and quantum trajectories 
solutions is taken from Ref. [25] (A = V/2, V = 57 ). 


B. Infinite dimension limit 

The qualitative differences between the variational ap¬ 
proach and the hierarchy equation method warrants fur¬ 
ther discussion, especially regarding the appearance of 
the bistable region. Some previous works have inter¬ 
preted this region as a genuine thermodynamic phase 
In this context, the concept of the lower criti¬ 
cal dimension is particularly important. It refers to the 
spatial dimension of a system above which a phase tran¬ 
sition can be observed [36]. For example, the equilib¬ 
rium liquid-gas transition belongs to the Ising universal¬ 
ity class and has a lower critical dimension of one. 

In the following, we will investigate the dissipative 
Rydberg gas in the limit where the coordination num¬ 
ber z goes to infinity. In this case both the variational 
method and the first-order hierarchy equations become 


exact, as the true steady state of the master equation 
will be given by a product state. In particular, it is 
instructive to look at this limit to investigate the role 
of the bistability found in the solution of the hierarchy 
equations. For this, we analyze the residual dissipation 
between the two local minima found by the variational 
method. As shown in Fig. [5] the variational method al¬ 
ways yields a unique steady state. Here we find that the 
larger the coordination number becomes, the smaller the 
difference in residual dissipation. However, the regime 
of true bistability indicated by a vanishing of the slope 
of the variational norm is only reached asymptotically as 
1/z, and for any finite value of z (i.e., for finite spatial 
dimensionality), there is no bistability. Consequently, it 
is incorrect to interpret the bistable behavior predicted 
by the hierarchy equations as a signature of a genuine 
thermodynamic phase, as it does not have a finite lower 
critical dimension. 
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FIG. 4: Extension of the antiferromagnetic phase. The 
shaded area depicts the presence of antiferromagnetic order 
according to the variational approach (left) and the first-order 
hierarchy equations (right). (V = 57 ) 


Rather, these findings suggest that the variational 
method is the correct starting point from which argu¬ 
ments in favor of the existence of thermodynamic phases 
in sufficiently high dimensions can be based. This is of 
course in contrast with equilibrium systems, where such 
arguments can be made based on a mean-field decou¬ 
pling, which is the equivalent of the first-order hierarchy 
equations. 


It is worth mentioning that the situation is very dif¬ 
ferent when the hierarchy equations predict an antifer¬ 
romagnetic phase. There, we also find an antiferromag¬ 
netic phase within the variational method, albeit with a 
smaller extension. While these results puts the existence 
of such an ordered phase on firmer grounds, the role of 
quantum fluctuation could still preclude its observation 
in actual experiments, if the lower critical dimension of 
the transition is three or larger. 
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FIG. 5: Analysis of the dissipative Rydberg gas in the limit of large coordination number z. The left panel shows the position 
of the first-order transition indicated by the crossing of the two local minima in the variational norm. The right panel displays 
the slope of the variational norm right below the transition point. The solid line is an algebraic fit to the data with an exponent 
of —0.99 ± 0.01, perfectly consistent with a 1/z behavior. 


C. Superatom model 


per superatom to be 


We now extend the previous discussion of the lattice 
model of Sec. IV Al to the case of a continuum, as it has 
been the case in the experimental situation in Ref. IH- 
In such a case, the Rydberg blockade will ensure that 
Rydberg excitations spontaneously form ordered struc¬ 
tures [33, Hij]. Although these correlations are short- 
ranged, we may still well replace the underlying contin¬ 
uum by a lattice structure with a lattice spacing corre¬ 
sponding to the typical spacing between Rydberg exci¬ 
tations. We can then determine the lattice spacing in a 
self-consistent manner [39| . finding 

= V^e ff + (2A) 2 . (17) 

The factor of two in front of the detuning A can be under¬ 
stood as realizing the antiblockade condition Ce/a 6 = 2A 
in the limit of vanishing Rabi frequency. The effective 
Rabi frequency fl e ff is derived from a renormalization of 
the atomic Rabi frequency 12 due to a (limited) collective 
enhancement. Far away from resonance, the transition to 
the first Rydberg excitation is still collectively enhanced, 
but the second Rydberg excitation can then only appear 
at specific positions that satisfy the antiblockade con¬ 
dition. Assuming there is always exactly one distance 
where the antiblockade condition is fulfilled, we find that 
we can describe the dynamics of such a superatom in 
terms of the number of atoms inside the superatom, N s , 
by renormalizing the Rabi frequency as 12 e g- = 
i.e., the geometric mean of the Rabi frequency for the 
first and the second Rydberg excitation. 

In the following, we assume the underlying lattice to 
be a cubic lattice (i.e, z = 6), however, we would like to 
stress that our results are basically independent of z, as 
long as it chosen consistently with the assumed lattice 
structure. Then, we can compute the number of atoms 


N„ = ‘ 




2VA 2 +12 2 ’ 


(18) 


with n being the density of ground state atoms. In the 
experimental situation of Ref. fl9| . the system was ei¬ 
ther on resonance or far away from it, i.e., either the 
condition fl A or 12 <C A has been fulfilled. Addition¬ 
ally, we capture the experimental situation by including 
a dephasing term associated with the laser linewidth n, 
according to the jump operators c' = y/2 7 nPr' > (40]. 

We are now in the position to compute the steady state 
of the system using the variational approach for corre¬ 
lated states. Remarkably, our results are in good quan- 
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FIG. 6 : Comparison of the stationary state of a dissipative 
Rydberg gas obtained by the Pisa experiment (left, taken with 
permission from 0 ) and by the variational method (right). 
For the experimental data, the color coding refers to the total 
number of Rydberg excitations in the sample, while for the 
numerical simulations, it represents the fraction of excited 
superatoms (7 = 5kHz, k = 500kHz, n = 1.8 x 10 17 cm -3 , 
C 6 = 7.54 x 10 - 58 Jm -6 ). 
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titative agreement with the experimental observations in 
Ref. [Il|, see Fig. [HI The jump in the density of Rydberg 
excitations corresponds to the first-order phase transition 
between a low-density gas and a high-density liquid, see 
Sec. IV Al In close analogy to the classical liquid-gas tran¬ 
sitions, the first-order transition ends in a critical point. 
However, due to experimental limitations on the dephas¬ 
ing rate and on the Rabi frequency, it appears that the 
observation of the critical region of the dissipative Ryd¬ 
berg gas remains a significant challenge. 

VI. SUMMARY 

In summary, we have given are detailed discussion of 
the recently introduced variational principle for st eady 
states of dissipative quantum nrany-body systems |25| . 
We have exemplified its usefulness by focusing on the 


driven-dissipative Rydberg gas, and we have made a com¬ 
parison of the variational approach to hierarchy equation 
methods, finding severe conceptual advantages in favor of 
the variational approach. Finally, we have found remark¬ 
able quantitative agreement with experimental data for 
the phase transition between a low-density gas of Ryd¬ 
berg excitations and a high-density liquid. Our results 
strengthen the position of the variational method as a 
key tool to analyze dissipative quantum many-body sys¬ 
tems. 
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